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Abstract 

Scattering power T = d < 0^ > /dx, used in proton transport calculations, is properly viewed as 
a differential description of the Gaussian approximation to multiple Coulomb scattering theories 
such as that of Moliere. Accurate formulas for T must take into account the competition between 
the Gaussian core and the single scattering tail of the angular distribution, which affects the rate 
of change of the Gaussian width, and leads to the single scattering correction. Mathematically, 
that implies that T must be nonlocal. In addition to proton energy and properties of the 
scattering material at the point of interest, it must depend in some way on how far the multiple 
scattering process has proceeded. 

We review five previous formulas for T and propose a sixth, 'differential Moliere', formula 
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TdM = fdMipV,PlVi) X — 

\pv J As 

where 

/dM = 0.5244 + 0.1975 lg(l - (pi; /piui)^) + 0.2320 Ig(pu) - 0.0098 lg(pt;) lg(l - (pv/pivif) 

is a fit to the single scattering correction as deduced from Moliere/Fano/Hanson theory. The 
scattering length Xs is a new material property similar to radiation length, pv is the product 
of proton momentum and speed at the point of interest, pivi the same at the initial energy, and 
Es = 15.0 Mev. TdM is easily computed and generalizes readily to mixed slabs because /dM is 
not explicitly material dependent. 

Whether or not an accurate formula for T is required depends very much on the problem 
at hand. For beam spreading in water, five of the six formulas for T give almost identical 
results, suggesting that patient dose calculations are insensitive to T. That is not as true of 
beam spreading in Pb. Evidently some favorable cancelation occurs in low-Z materials. At the 
opposite extreme, the projected rms beam width at the end of a Pb/Lexan/air stack, analogous 
to the upstream modulator in a passive beam spreading system, is sensitive to the choice of T. 
A simple experiment would discriminate between all but two of the six formulas. 

The usefulness of scattering power as a concept applies just as much to Monte Carlo as to 
deterministic transport calculations. For instance, using T in any of its forms will avoid step 
size dependence. Using the best available T might be important in general purpose Monte Carlo 
codes, which are expected to give the correct answer to many different problems. 



* Harvard University Laboratory for Particle Physics and Cosmology, 18 Hammond St., Cambridge, MA 01238, 
USA, bgottsch@fas.harvard.edu 
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1 Introduction 



Consider an initially monoenergetic, monodirectional proton beam slowing down in a finite slab of 
some homogeneous material. At any depth x we can write down the stopping power S = —dE/dx, 
the rate of decrease of energy with depth. There exists a well-known theory [1] for computing S, 
which depends only on the mean proton energy at x and atomic properties of the stopping medium. 
Integrating S over x we obtain, accurately, the total energy loss in a finite slab. 
At any x we can also calculate a scattering power 

T = d < 9^ > /dx (1) 

the rate of increase with x of the mean squared projected multiple Coulomb scattering (MCS) angle0 
There are a number of formulas for T in the literature. According to some of these T, too, depends 
only on local variables at x (mean proton energy and atomic properties of the stopping medium). 
However if we integrate these over x we do not, in general, obtain the correct rms MCS angle. The 
analogy with S is flawed. 

This might seem irrelevant since there exist accurate and well tested theories, notably that 
of Moliere [H 131 which directly give the total MCS angle, at arbitrary incident energy, for 
homogeneous slabs of any element, compound or mixture from very thin to near stopping thickness. 
Moreover, the derivation [51 [3] of Moliere theory does not rely on the concept 'stopping power'. It 
is never mentioned. 

Seen in that light, the aim of this paper is to discuss differential approximations to Moliere 
theory, functions T which, if integrated over x for any slab, will recover the Moliere MCS angle more 
or less accurately. 

Why do we want to do that? Because practical problems in proton transport go far beyond 
merely finding the MCS angle in a finite slab. We might for instance wish to compute transverse 
spreading of a pencil beam in a single slab, or in a stack of slabs, or the equivalent source poinlH 
in a stack, or the transverse penumbra of a proton therapy beam formed by such a stack. Such 
computations are done either by Monte Carlo methods or by deterministic (numerical) methods 
using the Gaussian approximation to proton transport, generalized Fermi-Eyges theorylf] Either 
way, we need to know, at every depth, the rate of change of the mean squared MCS angle. That 
need is exemplified by the formulas 

Anix) = I {x - x')''T{x')dx' , n = 0,l,2 (2) 

for the Fermi-Eyges moments, but it is more fundamental than that. Any transport calculation 
ultimately needs a recipe for the local rate of change of whatever variables are involved. 

It should go without saying that, when we do integrals like ([2]) numerically, approximating them 
by sums as in the midpoint rule or Simpson's rule, the answer should not depend on the step size Ax 
over some reasonable range. The same applies if Ax happens to be the step size in a Monte Carlo 
calculation. We mention this only because, in the case of T, this issue has lead to some confusion. 
For instance, Li and Rogers [8J in an electron papeiQ discuss the 'slab thickness dependence' of T 
at length, while a proton paper by Russell et al. [9] uses, without much discussion, Moliere theory 
with the step size (there conflated with slab thickness) fixed at an arbitrary value. 

^ In the early literature T is the rate of increase of the mean squared space angle, which is greater by a factor 2. 
In transport calculations projected quantities are usually more convenient. 

^ We use the term loosely here. Various 'equivalent sources' are defined in the literature cited below. 

^ We assume some familiarity with generalized Fermi-Eyges theory. For reviews see, for instance, ICRU Report 35 
[5] or more recent papers by HoUmark et al. [6] and Kanematsu [7_. 

* For historical reasons Fermi-Eyges theory and the concept of T were first developed (in the radiotherapy context) 
for electrons, but the conceptual issues are the same for protons. Indeed, protons are much better Fermi-Eyges 
particles. The small angle approximation comes naturally and the correlation of mean energy with depth holds to far 
greater depth because range straggling and the MCS detour factor are much smaller. 
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We shall see that there is, in fact, no problem provided we keep the concepts 'step size', 'slab 
thickness' and 'depth of the point of interest (POI)' carefully separated. All numerical results given 
later hold over a wide range of step size. We often use Ax (cm) equivalent to 0.1 g/cm^ but a much 
smaller or larger step gives essentially the same answer. 

The present paper deals with T and its implications in mixed slab geometry, with the slabs 
assumed to be infinite transversely. In most practical problems, whether beam line design or dose 
distribution in the patient, transverse heterogeneity or beam limiting devices eventually come into 
play. We must go beyond Fermi-Eyges theory, which is then relegated to propagating the individual 
pencil beams which make up the final (non Gaussian) distributions. The present paper ignores all 
that. We treat only the essentially one dimensional mixed slab problem, in order to understand 
that fully before introducing complications. Nor do we consider non electromagnetic effects like 
the halo from non-elastic nuclear reaction secondaries [TT] . This is strictly about the Gaussian 
approximation to MCS in mixed slab geometry. 

This paper overlaps a recent one by Kanematsu , 7: , whose notation we use except as noted. That 
work discusses the entire transport problem in the framework of Fermi-Eyges theory, treats heavy 
ions as well as protons, introduces a new scattering power TdH and presents many useful analytic 
approximations. The present paper focuses entirely on scattering power, protons, and single or 
mixed slabs, largely glossing over computational issues. It is an in-depth look at a small part of 
the ground covered by Kanematsu. His work is somewhat biased towards dose reconstruction in the 
patient, ours towards beam line design, but the fundamental issues are the same and there is no 
significant disagreement. 

For completeness, we repeat a number of formulas from the literature, some of which is quite old. 
There is, however, considerable new content: simplifying the formula we call Tjc, deriving a new 
non-local formula T^Mj and proposing, for a case of practical importance, an experimental test to 
discriminate between various formulas for T. In short, we combine original material with a critical 
review which, we hope, will be useful to the student and of some interest to the expert. 

2 Preliminaries 

2.1 Mixed-Slab Notation 

Take, for concreteness, an example from beam line design: a monoenergetic proton beam of known 
emittanc^l enters a Pb/Lexan/air stack (Figure [1]). That is precisely what happens at the upstream 
range modulator/first scatterer of a modern passive beam spreading system [12 |f| (In the patient, 
the variation of material properties is usually less extreme, but still present.) We wish to transport 
the beam through the stack to an arbitrary measuring plane MP at a depth x (cm) measured from 
the entrance of the first slab. In particular let us assume we wish to find yrms of the Gaussian which 
approximately describes the transverse fiuence distribution on an MP located, say, at the second 
scatterer, because the precise match of yrms to the second scatterer design [13 will determine the 
fiatness of the dose distribution at the patient. 

Mi is the material of which the i**^ slab is composed and stands for a set of material properties 
e.g. density pi (g/cm^), radiation length Xoi (cm) and (to be defined) scattering length Xsi (cm). 
There are also atomic weights A, atomic numbers Z and fractions by weight w of the constituents of 
Mi, all implicit stepwise constant functions of x. Xi refers to the entrance of the i^^ slab and xi = 0. 

Ei (MeV) is the proton kinetic energy entering the i*^ slab, piVi (MeV) the corresponding product 
of momentum and speed, and Ri (cm) the mean proton range corresponding to Ei stopping in M^, 
that is to say the residual range in Mj at the entrance to the i**^ slab. E, pv and R (the latter always 
computed for the current material) are corresponding quantities at the depth of interest x. Given a 

^ In other words the parameters of the incident Fermi-Eyges beam eUipse [5] are known to sufficient accuracy. 
" An 'upstream' modulator also serves as the first scatterer in a double scattering system, and MCS in it is critical. 
A 'downstream' modulator is near the patient and MCS has little effect. 
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specified stack and any value of x witliin it, all these quantities can be computed using range-energy 
tables and kinematic relations^ 

This example is a rare practical case where Fermi-Eyges theory alone gives a useful answer, 
because there is in fact no transverse heterogeneity and the final fluence distribution is in fact 
very nearly Gaussian. In more complicated beam line design problems, with collimators as well as 
scatterers whose thickness varies with radius, Fermi-Eyges theory can only be a building block, as 
noted earlier. 

2.2 Limits of Moliere Theory 

Radiotherapy protons have kinetic energies roughly in the range 3 < E < 300 MeV corresponding to 
0.015 < i?i < 51.5 cm in water. The interesting range of normalized depth in a single slab is roughly 
0.001 < x/Ri < 0.97. The lower limit might for instance apply to a vacuum window in a beam 
scanning system. The upper, near stopping depth, is where Moliere theory fails 4J as straggling 
blurs the correlation between depth and energy. It is reached and exceeded when protons stop in 
the patient. 

It is well known that Moliere theory applies only to slabs of some minimum thickness. There 
must be enough atomic encounters to establish the regime of multiple (as distinct from single or 
plural) scattering!! Moliere .2j gives B — 4.5 as the lowest allowed value of his slab thickness 
parameter Ifl Table[T]lists target parameters at which that limit is reached for some common materials 





P 


pRi 


X 


x/Ri 


^Hanson 




g/cm"^ 


g/cm^ 


/im 


xlO^ 


mrad 


Be 


1.85 


21.29 


0.65 


5.6 


0.025 


water 


1.00 


17.38 


2.35 


13.5 


0.051 


air 


0.0012 


19.67 


0.24 cm 


14.4 


0.056 


Cu 


8.96 


26.26 


0.79 


27.0 


0.167 


Pb 


11.35 


36.06 


2.11 


66.6 


0.483 



Table 1: Target parameters for which Moliere's B equals 4.5 at 158.6 MeV incident. 

for 158.6 MeV incident protons. It is clear that the lower Moliere limit need not concern us in 
practical proton radiotherapy calculations. 

2.3 Kinematics 

The kinematic expression 1/pv appears in all multiple scattering formulas If we define the reduced 
kinetic energy of any particle as 

E 



(3) 



where mc^ is the particle's rest energy, we find 



pv = ^E (4) 

T + 1 



^ In the early literature x is frequently expressed in units of the radiation length Xq (therefore dimensionless) and 
instead of T, the mass stopping power T/p (radian^ /(g/cm'^)) is used. This is unhelpful for mixed slabs. 

* One can speak of the (extremely small) average energy loss in an atomic monolayer, or of the (extremely small) 
single scattering probability, but it makes no sense at all to speak of multiple scattering since there is at most one 
collision. That is the fundamental difference between stopping and multiple scattering theories. 

^ B (X logarithm of normalized slab thickness. The constant of proportionality varies with the material [4]. 

Multiple scattering is derived from single scattering as outlined in Scction FS.ll 1/v enters the derivation of the 
single scattering probability I I14I I because the impulse delivered in a single collision is proportional to the interaction 
time or inverse speed. 1/p comes from the fact that the angle of deflection equals Ap/p. 
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For 3 < E < 300 MeV protons the coefficient of E varies from 2.00 to 1.76. Other useful relations 
in the same vein are 

and 

(pc)2 = {T + 2)mc^E (6) 

These formulas avoid differences between large quantities that can arise in relativistic calculations, 
and their relativistic and non-relativistic limits are obvious by inspection. 

2.4 The 0veras Approximation 

0veras |14| found a simple empirical relation between {pv)'^ and normalized residual range. For a 
single slab 

(H' = (Pi^'i)' (7) 

where k depends on the material. If the material is characterized by its radiation length, Schneider 
et al. jl5] found the empirical expression 

fc = 0.12 e"^-^^'^^" + 0.0753 (8) 
k is fairly small compared to 1, and often the 'weak 0veras' approximation 

(H' ~ (Pivi? a-x/Ri) (9) 

is adequate. We shall use it later as a guideline. Figure [2] tests ^ over a wide range of materials, 
normalized slab thicknesses and outgoing kinetic energies. 

2.5 Moliere/Fano/Hanson Procedure: 6'Hanson5 ^Hanson 

For a single slab, given the outgoing energy, material and thickness, we can find the rms angle in 
the Gaussian approximation by using the Moliere/Fano/Hanson procedure [4] with 

'^Hanson 

{E,Mi,x) = XcV'B - I.2/V2 (10) 

We use this angle as our standard of comparison for all rms angles obtained by integrating scattering 
powers over single slabs. It is Kanematsu's ^mh jZ]- We have given it the longer name to emphasize 
that it is not derived from a scattering power. Indeed, we now derive a scattering power from it by 
numerical differentiatiorF^ 

^^Hanson _ ^ ( j, ^ ^Ln.on(^. M,,x)- ^Lnson(-g. M^,X- Ax) 

J = J^HnnsoniE, Ml,X) = llm (iij 

ax Ax^o Ax 

Though we call this Tfjanson it is not a scattering power in the usual sense. Evaluating it requires a 
lengthy procedure rather than a simple formula. However it is the 'correct' T for single slabs and 
we will use it later to derive an improved scattering power 



We use the single sided derivative because adding even a small increment to x gets us into trouble for near 
stopping length slabs, where ^Hanson levels off |4] and Tjjanson = 0. 

Eq. mil l may be compared with Eq. (11) of Li and Rogers [S] which, it seems to us, defines an average rather 
than an instantaneous rate of change. 
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2.6 The Single Scattering Correction 



Rather than thinking of a proton entering a given slab at some energy, let us focus on the proton 
at a depth x or equivalently, a proton leaving a single slab of thickness x. Assume the material is 
Be and the energy at x is 20 MeV. A 'local' formula for T uses only those two facts. But there is 
another parameter: the amount of material overlying x. Our 20 Mev might for instance represent an 
incident 23.7 MeV proton passing through an overlying 1mm of Be {x/Ri = 0.264), or a 102 MeV 
proton passing through an overlying 5 cm {x/Ri = 0.95). Obviously the mean squared MCS angle 
itself depends not only on E{x) and M{x) but on the quantity of 'MCS buildup' material x/Ri 
traversed to get to x. Does T, the rate of change of mean squared angle with a;, also depend on 
x/Ri? 

It does. Figure [3] shows mass scattering power T/p, computed three different ways for three 
materials at 20 MeV, as a function of normalized overlying material Tkanson IHl) is the 

'correct' scattering power. It is obviously nonlocal: it depends on x/Ri. Tfr and Tic (discussed 
later) are local: they do not depend on a;/i?i. Of the two, Tic has the more accurate material 
dependence suggesting that it, with some logarithmic function of x/Ri for the single scattering 
correction, might lead to an accurate approximation to Tfjanson- 



2.7 Highland's Formula: ^Highland 

Highland |16j . in order to simplify MCS calculations for the high-energy physicist, parameterized 
Moliere/Bethe/Hanson theory [3] and obtained the elegantly simple formulcO 



14.1 MeV f 1. X 



logio ^ (12) 



in which the single scattering correction, in parentheses, is evident. He assumed a slab that is finite 
but sufhciently thin that pv does not decrease much (x <C Ri). Thus p2|) is already an integral 
expression. It cannot be applied to a thicker slab by interpreting a; as a step size Aa:, using (fT^ for 
each step, and adding in quadrature: the sum decreases indefinitely with the number of steps. 

We circumvented that problem in [1] by arbitrarily removing the logarithmic term from the sum 
(or integral) and proposing a generalized Highland formula 

^Highland = (l + ^logio ) I / ( ^ ) dx' ] (13) 

which reduces to (fT^ for thinner slabs and was found [3] to agree with measurement almost as well 
as Moliere/Fano theory. It is Kanematsu's 9iH , which we have given a longer name to emphasize that 
it is not computed from a scattering power. Figure[3] compares 6'Highiand (fT5|) with ^Hanson (fTU]) . The 
behavior for Be reflects the difference between the Moliere/Bethe/Hanson and Moliere/Fano/Hanson 
forms of the theory. Otherwise, ^Highland is better than the ±5% advertised |l1^, and much easier 
to compute than the full theory. 




3 Prior Formulas for T 

At this writing, three local and two nonlocal formulas for T can be found in the literature. We 
review the first two at some length because Rossi's excellent book is no longer easily available 
and because the second. Tic, is the basis for our improved T. 

See Appendix|X]for details on how we computed examples, figures and tables. 
1"* The 'Eb constant', strictly following Highland's paper, should be 17.5 X 1.125/^2 = 13.92. However, the 1986 
Particle Properties Data Book gave 14.1 which, for whatever reason, we [1] used in our comparison with experimental 
data. That has since become the accepted value, used by both Schneider et al. [15] and Kanematsu [?]• 
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3.1 Fermi-Rossi: Tfr, 6'fr 

Following Rossi Q7] the single scattering probability for a singly charged particle in small angle 
approximation is l£I 

S(x) dndx = 4 Nrl ^ ic^) ^ ^dndx (14) 

Te is the classical electron radius, N is Avogadro's number, A, Z and p are the atomic weight, atomic 
number and density of the target material and rrieC^ is the electron rest energy. 

In the derivation of (HH), the target nucleus is modeled as an unscreened point charge. More 
realistically, the scattering law must depart from 1/x^ at small angles (distant collisions) of order 

XI = 1.13aZi/3 (^a^^ (15) 

due to screening of the nuclear charge by atomic electrons, and at large angles (close collisions) of 
order ^ 

= ^ (^) 

due to the finite size of the nucleus, a is the fine structure constant. 

Rossi now assumes that the value of 0^ at a; + dx equals its value at x plus the mean squared 
space angle of scattering in dx. This step is equivalent to assuming the process is exactly Gaussian 
and ignoring the single scattering correction. It leads to 

de^ = pdx / x' s(x) dn^ (17) 

Jo Jxi 

Rossi now defines 

^l^z^ - / x'^ix) dn^ (18) 



p dx 



XI 



Later, Brahme 18 used [T/ p) and the term 'mass scattering power' for 0^. We remind the reader 
that this and other early uses of (T/ p) refer to the rate of change of space angle unlike later treatments 
and the present paper. 

To perform the integral in (|18|) analytically we must assume some simple behavior of S(x) below 
Xl and above X2- Rossi does this two different ways. In the first, Rossi assumes S is zero in both 
regions. After integrating, simplifying and introducing^ 



27r 



1/2 



Es = y— ) ruec^ = 15.0 MeV (19) 
he obtains 

T^FR = f^V 4;^ (20) 



pv J X, 



Xq (cm) is the radiation length of the material. For a sufficiently thin single slab pv does not change 
much, the integral 

Ao =<e^ >= I T{x')dx' 







^® We use Moliere's notation x foi' single scattering space angle and related parameters and, in order to follow Rossi 
more closely, briefly use for the multiple scattering space angle to distinguish it from 9, the projected angle. 

^® We have kept two changes from ICRU Report 35. The 1.13 comes from the Thomas-Fermi radius of the atom, 
Ta = 0.885 ao Z"^/^ (ao = Bohr radius), where Rossi used 1 instead of 0.885. In Eq.JTOj ICRU35 rounded Rossi's 
(1/0.49) to 2. 

Here we switch back to projected angle 8. 
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is trivial and 

^™ - (21) 

the well known Rossi formula. For thicker slabs, pv can be related to x' as long as we know the 
range-energy relation in the material, and the integral is performed numericallvF^ Figure [5] compares 
^FR, so computed, to ^Hanson- ^FR IS far too large for thin scatterers, by an amount which depends 
on material. 



3.2 ICRU Report 35: Tic, ^ic 

Rossi next cites a much better approximation, assuming that 5(x) behaves as l/(x^ + Xi)^ small 
angles (leveling off at small x rather than suddenly vanishing) . Eq. ^TEh then gives 

^ aNrli^X ^{\og[l+[^] I - 1 + I 1 + ( ) I V (22) 





pv J A 

This is in essence the scattering power given in ICRU Report 35 [5], except that the version given 
there takes advantage of a chance cancelation with WeC^ and therefore, unlike ()22|) . only applies to 
electrons. Nevertheless, we shall call it Tic to distinguish it from TpR. 

For protons. Tic can be simplified considerably. When X2 (|16p comes out larger than 1 radian it 
should be truncated to 1 [17]. For radiotherapy protons that never happens. The worst case is 3 Mev 
protons in Be where X2 = 0.913 rad. We can therefore simplify X2/X1J canceling the p dependence. 
If we also ignore the rightmost term in X2/X1J which is always much less than 1, and introduce a 
scattering length defined by 



1 . o Z 



aNr 



|21og(33219 (AZ)-i/3) - l| (23) 



we find for radiotherapy protons, 3 to 300 MeV, 

Tic i- (24) 

\pvj Xs 

identical in form to TpR (pO|) 

For compounds or mixtures, any scattering power obeys a Bragg rule. Atoms act independently, 
and the compound or mixture is equivalent to very thin sheets of each constituent in the correct 



compi 

proportion!^ That picture leads to 



pXs ^ \pXs 



where Wi is the fraction by weight of the i"^ constituent. Xq obeys a similar formula. Table [2] 
compares Xs with Xq for a few materials. Figure [5] compares ^ic, obtained by integrating Tic, with 
^Hanson- Material dependence is greatly improved over TpR, but the large error for thin scatterers 
remains. 

HoUmark et al. [B] use Tic for protons and heavy ions. However, the formula they quote (their 
Eq. (24)) is only valid for electrons, and they introduce an effective charge factor Zp^cS without 
commentl^ Their values of T/ p for protons up to 200 MeV in water (their Table 2) are larger than 
ours by a factor 1.19, perhaps due to incorrect application of the Bragg rule. Our Figure [T^ seems 
to confirm that discrepancy. A footnote to a more recent paper [llj by the same group corrects the 
headings of Table 2 but does not mention any numerical error. 



See Appendix lAl for details. 

On fundamental grounds one would expect the Bragg rule for scattering powers to be very much better than for 
stopping powers, where (through /) there is some sensitivity to molecular binding [I]. 

^'^ If we use the Barkas formula Z^g = Z{1 — exp(— 125/3Z~^/^)) [T9l|20] then Zp^eff = 1 for radiotherapy protons. 
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pXs (g/cm^) 
pXq (g/cm^) 



Be Lexan 

92.60 55.05 

65.19 41.46 

1.420 1.328 



H2O Al 

46.88 28.75 

36.08 24.01 

1.299 1.197 



Cu Pb 

14.62 6.62 
12.86 6.37 

1.137 1.040 



Table 2: Comparison of scattering length Xs with radiation length Xq for six materials. 
3.3 Linear Displacement: Tld, 6'ld 

In a recent note [21j Kanematsu proposes a simple scattering power for protons and heavy ions in 
tissue-like matter. He uses water as a reference material. For protons, his formula reduces to 

Tld = 1.00 x 10"^ ^ (26) 

In a single slab of tissue-like material M 

Px = Xqyj/Xqm (27) 
Rw = Riw - PsX (28) 
Ps = Sm/Sw (29) 

where W stands for water, S is stopping power (MeV/cm) and Rw (cm) is the proton's residual 
range in water at depth x. In the interest of a more complete survey we will, in what follows, explore 
the behavior of Tld for non tissue-like materials, even though it may not work well for those. We 
call it the 'Linear Displacement' scattering power after Kanematsu's derivation. For a single slab it 
can be integrated analytically giving 



Ao ^<9l^>= I TMx')dx' = 1.00 X 10-3 J. In ) (30) 



^0 



A plot of VS. x/Ri (Figure [7]) is reminiscent of 0fr except near end-of-rangeF^ Indeed Tld is a 
variant of Tp^ with three changes: the weak 0veras approximation ^ is used for pv, Eg is adjusted 
downward, and water is used as a reference material. To show this consider a single slab of water. 
Then 



Eg 1 Eg Riw 1 



IFR.W 



frt 



^ow {pv^ Xow ipv)l Riw - X 
152 I I 

X (2 X 10^^) = 1.25 X IQ-^ — (31) 



36.08 ' ' Rw Rw 

if we evaluate -Riw/(pi'^i)^ at 158.6MeV (it is insensitive to Ei). When 1.25 is reduced to 1.00 and 
([3T|) is generalized to other materials by introducing ratios to water, (|26| follows. 

Kanematsu [21 freely admits that Tld is semi-empirical, but his derivation [7ll21j is complicated 
and it is somewhat unclear exactly where the downward adjustment of Eg takes place. It is probably 
mostly present in his Eq. (6) |21j which flows from Highland's equation which in turn is an empirical 
fit to Moliere/Bethe/Hanson theory 16 . 

There is precedent for adjusting Eg downward to improve the performance of TpR over a limited 
range of normalized thickness. Soukop et al. multiplied Eg by a factor 0.8, said [22] to have been 
obtained from a fit to Geant 4, in their 'corrected Rossi' formula. Kanematsu's reduction factor is 
only (1/1.25)1/2 ^ o.89. 



For the sake of uniformity, Figure [7] was obtained by integrating II26I I numerically. A small difficulty arises for 
non tissue-like material because the stopping power ratio (I29I I is then somewhat energy dependent. That introduces 
some irregularity near end-of-range, which we minimized by evaluating ps at 0.64 X Ei . 



9 



3.4 0veras-Schneider: Tqo,, 6*03 



Schneider et al. |15j describe a scattering power based on Tfr with a nonlocal correction factor in 
the form of an analytical function of a new normalized variable t. For a single slab, t[x) is just x/ Ri. 
For mixed slabs 

t{x) = x/R{Ei,M{x)) (32) 

where M{x) is the current material, i is a discontinuous function of x. It is that normalized depth 
which would be obtained if the protons were degraded from Ei to E{x) in the current material. 

They derive their function of t by fitting a large body of single slab experimental data for Orms 
at 158.6 Mev [4] with a two parameter analytic function of t (their Eq. (10)). Differentiating, they 
fincS 

/19.9MeV\^ 1 

V PlVi J Xq 

with fitted constants fc(Xo) ® and 

Co = (201/200) - (23/5000)pXo (34) 
Cl = -(11/2) + (43/1000)pXo (35) 

For mixed slabs, pX^ in Eqs. ([5]), ([M]) and ((55|) is the mass radiation length of the current 

material. Figure [5] shows considerable improvement over the local formulas for T (note change in 
scale). Oscillations and divergent behavior at the ends are characteristic of polynomial fits. 

Two general remarks. First, Schneider et al. fit experimental data rather than some form of 
multiple scattering theory. Since theory and measurement agree rather well [4] this should not have 
a large effect, but it biases their formula towards the data that happen to be available. Second, 
their major advance is the introduction of a nonlocal correction based on a generalized definition of 
normalized depth. The 0veras approximation, tightly woven into their formalism to obtain formulas 
in closed form, is somewhat of a distraction. Similar results could be obtained without it. 



3.5 Differential Highland: TdH, ^dH 

Kanematsu [7^ derives a nonlocal scattering power applicable to mixed slabs which by construction 
gives the same result as (IT^ for a single thin slab. In his case the nonlocality parameter is depth 
weighted by inverse radiation length: he generalizes x/ Xq to a dimensionless radiative path length 

r dx' 

i{x) ^ / (36) 
and writes a new scattering power as Tfr, times a correction factor 

TdH ^ /dH(^) (—)' ^ (37) 
\pv J Xq 

He now constructs / so that, for a single thin slab, the integral of TdH will equal ^Highland- That 
requires that the average of / equal the squared ratio of the Highland formula (fT2|) to the Rossi 
formula 

lg£\ /14.1MeV\^ 



We have corrected two typographic errors. 
23 Ig = login, In = log„. 
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Differentiating, he finds after some algebra 



fMi) = 0.970 (l + ^ j (l + ^ j (38) 

to be used in (jST]) . Figure ([9]) compares ^dH, obtained by integrating TdH, with (?Hanson- It indeed 
behaves very Hke ^Highland except for thick slabs because the derivation of ([55]) is strictly correct 
only for thin slabs. 

Unlike the extension of TdH to mixed slabs is totally straightforward, requiring only the new 
path integral £{x). 



4 Improved Nonlocal Formula: T^Mt Odu 

Instead of Tfr as a basis let us use Tic ([M)) which, for radiotherapy protons, is as simple and has 
better material dependence (Figure (HI). From Figure [31 the correction should be logarithmic in total 
material overlying the POL 

In the weak 0veras approximation ([9]) we found that 1 — (pv/pivi)'^ , which depends only on 
local energy E and incident energy ii^i, is a reasonably good proxy for normalized depth x/Ri for 
all materials and energies of interest over three orders of magnitude (Figure [2|) . Let us therefore 
compute and plot the ratio of Tnanson, the ideal scattering power, to Tic, our proposed basis, for 
0.001 < 1 - {pv/pivi)"^ < 0.97 and Ei < 300 Mev, for several values of the energy E at the POI 
and several materials (our usual Be, Cu and Pb). For each point we compute the exact x/Ri for 
that material, without relying on the 0veras approximation, then use (llip . The result is shown in 
Figure [ini 

An adequate fit to these data, also shown in Figure [TUl is a linear polynomial in Ig (1 — (pu/piui)^) 
whose two coefficients are in turn linear in Ig (pv). Defining a new 'differential Moliere' scattering 
power and writing everything out we have 

TdM = fdMipv,piVi) X ( — ) (39) 

\pv J As 

where 

/dM = 0.5244 + 0.1975 lg(l-(p«/pi«i)2)+ 0.2320 Ig(pv)- 0.0098 lg(pw)lg(l-(pw/piwi)2) (40) 

Although the weak 0veras approximation suggested the form of (|40p the final result is an independent 
fit. It does not depend on the accuracy of either ([7]) or ([§]), or on the scattering material. 

FigurefTT] compares O^m obtained by integrating TdM with 6'Hanson (note the vertical scale) . 

The reader will object that the coefficients of ([3D]), here given to excessive precision, are arbitrary. 
We could have chosen different sets of materials or energies. That is perfectly true, but it is also 
true of T0S and TdH- The former is a fit to a specific data set for specific materials at a specific 
energy. The latter stems from Highland's formula, also a fit albeit to theory. That said, we have 
tried other materials and energies without much change in the general appearance of Figure [TTl 

Since 0dM agrees well with ^Hanson it may be inferred [3] that it agrees well with experiment 
for many materials. We will show that directly only for one \ow-Z and one high-Z material, also 
using the opportunity for a head-to-head comparison of six scattering powers and the generalized 
Highland formula. Figure [T2l shows the comparison for polystyrene and Figure [T3l shows it for lead, 
both with data from [4]. All three nonlocal T's are better than any local T but TdM is the best. It 
almost agrees with measurement within the experimental error 

Reference [4] already remarked on the fact that experimental data seem a few percent lower than theory for thick 
Pb slabs. In our present opinion that is more likely due to a systematic experimental error for very large Brms than 
to a breakdown of Moliere theory, but it is impossible to say for sure. 
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To facilitate numerical checks we include short tables of 6'xx (Table S]) and (T/p)xx (Table [5]) for 
single slabs of various materials and normalized thicknesses. We use Ti = 158.6 MeV to correspond 
to particularly Table 1. ^Hanson is not given directly there but may be found using ([TI])) . 

The generalization of T^m to mixed slabs is the easiest of all. It does not even require an 
additional path integral. The single scattering correction is a logarithmic function of the fractional 
decrease in pv, with no explicit material dependence, from the incident beam to the point of interest. 
It diverges, as any single scattering correction must, at pv/{piVi) — 1 (no overlying material). 

5 Applications 

What kinds of computation are significantly affected by the choice of T? If we are only interested 
in the Gaussian MCS angle itself in a single slab. Figures [4l through [9l and fTD through [T3l alreadv 
answer that question. However, in that case we would not use transport theory or T at all but 
simply evaluate either 6'Hanson (UHl) or 6'Highiand (US])- Let us examine some less trivial cases. 

5.1 Pencil Beam Spreading in a Single Slab 

The archetype for dose reconstruction in the patient is pencil beam spreading in a homogeneous water 
slab. An early paper by Preston and Koehler "23J derived a universal formula for beam spreading 
and compared theory with measurements. HoUmark et al. [6] refer to additional measurements. 
Figure 5a of [4] already suggests that beam spreading in water is insensitive to T, comparing two 
very different models of MCS: Preston and Koehler's which is a local model similar to Tic, and a 
beam spreading model based on the nonlocal generalized Highland formula. Insensitivity of beam 
spreading in water to T implies that almost any T will work reasonably well for dose reconstruction, 
and is therefore worth a closer look. 

Figure [m shows spreading of a 127 MeV pencil beam in water with ymis = (^2)^^^ (Eq.[2|) 
computed according to all six formulas for T and according to Table 2 of HoUmark et al. [6j (Tho)- 
Except Tpji and ThOi all agree with experiment and are barely distinguishable from each other. 
Figure [15] is an expanded version where we plot the difference between each calculation and Preston 
and Koehler's formula (Appendix [B]) . The spread, of order 0.1mm, would be negligible in any dose 
reconstruction problem. 

To gain some insight into this insensitivity to T, Figure [TBI shows the integrand of A2 when the 
POI is at 0.97 near stopping depth. The near linearity of curves for the local formulas (Tfr, 
Tic and Tld) is a consequence of the 0veras approximation. The nonlocal formulas do in fact give 
a lower result because of the single scattering correction, important for small x, but not by much. 
Except for Tfr, which is too high everywhere, the areas under the curves are nearly the same. That 
breaks down when we consider the integrand to x = 0.1i?i (Figure fT7|) but by then the absolute 
effect is so small as to be negligible. 

It is instructive to look at the same problem in terms of the evolution of the beam phase space 
ellipse [5]. Figure [TSl studies the same case as the preceding figures, dividing the water slab up to 
X = 0.97 i?i into five sub-slabs. We also show bounding boxes for the final slabl^ The spread of 
vertical bounds shows that the different T's do give different answers for the final rms angle, but 
the effect on spatial spreading (horizontal bounds) is almost nil. The same study for Pb, Figure fT9l 
shows that this is a fortuitous property of tissue- like materials, presumably due to a particular 
combination of drift and scatter. Beam spreading in near-stopping high-Z slabs is not entirely 
academic. In computing collimator scatter, we are basically asking how many protons leave the bore 
of a Cerrobend or other aperture that, except for MCS, would have stopped. 

^® Recall that the vertical bound represents 9rms while the horizontal bound is yrms [5]. 



12 



5.2 MCS Angle in a Double Slab (Range Modulator) 



Unlike beam spreading in water, the choice of T is important in at least one practical problem in 
proton transport. The upstream range modulator in a passive beam spreading system is a sequence 
of high-Z/low-Z sandwiches designed to produce a Gaussian fluence distribution of constant width 
Vrms either at the patient or, in a double scattering system, at the second scatterer [T^HS]. Also, 
each sandwich is designed to pull back the pristine depth-dose by some fixed water equivalent amount 
so as to produce the desired spread-out Bragg peak. 

The design procedure [H] is fairly complicated and the details need not concern us. We simply 
wish to define a sequence of Pb/Lexan/air stacks that might occur in a practical beam line, transport 
the beam through each such stack, and see how much difference T makes. We could close the loop 
with an experiment to see which T is best. 

Table[6]hsts the parameters of a typical design, feanson, 3 is the angle obtained by combining the 
Moliere/Fano/Hanson angles for Pb and Lexan in quadrature. Scattering in air is ignored. More 
important, the single scattering correction for Lexan is wrong because its MCS angle is calculated 
for a beam of energy E2 entering the Lexan de novo, ignoring the fact that some MCS has already 
taken place. For this and other reasons we do not expect the design and transport calculations to 
agree exactly for any choice of T. 

FigureHnishows Fermi-Eyges transport results for each Pb/Lexan/air stack, with an ideal incident 
beamo We used Kanematsu's finite increment form ([^ Eqs. 19-21), equivalent to integrating by 
the midpoint rule. The finite depth of the midpoint of the first step sidesteps the divergence of non- 
local T's at a; = 0. We discretize each slab separately with a single minimum step size parameter 
pAx = 0.1 g/cm^ for all materials, yielding typically tens of steps in Pb, a few hundred in Lexan, 
and one step for air. Increasing pAa; even by a factor 20 changes yrms less than a percent. 

We have not yet done an experiment to see what yrms these or similar setups actually produce, 
but Figure [20I is encouraging. The fact that our standard design procedure in fact produces a flat 
dose at isocenter (matches the second scatterer design) suggests that yrms — 3.5 cm is probably 
correct to a few percent. If, lacking a direct measurement of yrms, we assume that to be the case, 
we find that T^m indeed comes closest to the right answer whereas some of the other T's are off by 
an amount much larger than the 1-2% a careful experiment could measure. 

It is amusing that Tld also gives a very good result, even in a non tissue- like problem. That is 
consistent with Figure[7l Orms for a moderately thick Pb slab is slightly low, whereas for a thick 
plastic slab it would be slightly high. It augurs well for 'corrected Rossi' T's in general [22] even 
though the excellent result here may be somewhat accidental (or, depending on the outcome of a 
measurement, possibly wrong). 

One could argue that, because we have an adequate 'slabwisc-Hanson' design procedure, we 
again don't really need either beam transport or T. And indeed we don't, if we are merely designing 
beam spreading systems. However, the Pb/Lexan/air stack could be part of a larger Monte Carlo 
computation [25j . That brings us to the last topic. 

5.3 Monte Carlo Calculations 

Monte Carlo calculations are the gold standard in radiotherapy. A condensed history Monte Carlo 
should embody a differential model of multiple scattering analogous to scattering power, even if it 
is not called that. In each (finite but small) step of known material, one needs to compute the 
increase in the width parameter of some distribution (Gaussian, Moliere or other) from which a 
random deflection angle is then drawn. Just as in deterministic calculations, one would like the 
final result to be independent of step size over some reasonable range. Table [3] shows for a simple 
case that results from Geant4, whose MCS model is based on a variant of Highland's formula [2^ . 
depend somewhat on step size. Monte Carlos based on a 'corrected Rossi' formula [2^, similar to 
Tld, should not have this problem, but may not always give the right answer (Figure[7]) since they 

A known initial beam phase ellipse could easily be used if necessary. 
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STEPMAX 


# steps 


^rms 


(^rms /^Hanson 






millirad 


% 


default 


12.1 


109.92 


-1.1 


O.lx t 


19.1 


111.13 


0.0 


O.Olx t 


101.6 


118.35 


6.5 


O.OOlxt 


1015.0 


123.19 


10.9 



Table 3: A test using Geant4 v9.1, courtesy L. Urban, CERN. 158.6 MeV protons scatter in 20.196 
g/cm^ Pb {Orms,cxpt = 108 ± Imrad, 6'Hanson = 111.1 mrad). STEPMAX governs the step size and 
Orms is obtained by fitting the projected angular distribution with a Gaussian. 

lack the single scattering correction. T^m, easily generalized to other charged particles, might be a 
good compromise between accuracy and step size independence. 

6 Summary 

Unlike stopping theory, accurate theories of multiple Coulomb scattering such as Moliere's do not flow 
from a differential form. If that is needed, for deterministic or Monte Carlo transport calculations, 
it must be devised retroactively as an approximation to the more exact theory. To be accurate for 
thin scatterers it must include a single scattering correction. That implies nonlocality: in addition 
to energy and properties of the material at the POI, it must depend in some way on how far the 
multiple scattering process has advanced. That does not conflict with being numerically integrable 
in the usual sense that an approximating sum approaches a limit as the step size decreases. 

Nonlocality may be characterized in different ways. Schneider et al. use a generalized definition 
of the normalized depth of the POI, Kanematsu uses a radiative pathlength integral up to the POI, 
and we use the diminution of pv from its incident value to the POI. All three can be tested against 
Moliere theory for uniform slabs, and all three generalize to mixed slabs, as they must to be useful. 

We have reviewed three local and two nonlocal formulas for T, comparing 9rms from each one 
graphically with the 'correct' MCS angle feanson- One of them, Tjc, can be simplified for protons 
by introducing a new material parameter, the scattering length Xs, which is similar in form to 
radiation length Xq. Tic provides a basis for a new nonlocal scattering power T^m (Eqs. l23l25|39l 
and HP]) which, for single slabs with 0.001 < x/Ri < 0.97, reproduces Sanson as well as the 'correct' 
scattering power Tnanson to a few percent. 

In practical problems, the choice of T is frequently not critical. We have shown, for instance, that 
all T's described here except Tfr work well for beam spreading in water and, presumably, water-like 
materials. In particular Tld, typical of 'corrected Rossi' formulas, works well. Beam spreading in 
Pb and, presumably, other high-Z materials, is more sensitive to T. 

Turning to a mixed slab problem, we have shown that Pb/Lexan/air combinations typical of an 
upstream range modulator yield very different answers for different T's, so that a direct experimental 
test should be easy. Again in this case, Tld (even though used outside its supposed range of validity) 
yields almost the same answer as T^m, probably because of a lucky cancelation. 

In the end, perhaps the strongest case for an accurate scattering power such as TdM can be made 
for general purpose Monte Carlo codes which are supposed to do everything well, rather than beam 
design or patient dose computations where special purpose workarounds can be devised. 
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A Computational Details 

Examples, figures and tables were computed with Fortran programs which may be downloaded free: 
see \BGware.zip at http://physics.harvard.edu/~gottschalk. Source code is in \BGware\source. All 
calculations are single precision. 

Formulas for variants of Moliere theory are given in implemented by module THETAO.FOR. 
The only significant change is that we now use cubic spline interpolation of lni?(lnTj^, rather than 
a polynomial fit, to interpolate range-energy tables (RANGE. FOR) Subroutine TOUT computes 
the energy out of a stack of slabs, as well as the outgoing rms projected MCS angle using what we 
have called the slabwise-Hanson or slabwise-Highland procedure, not Fermi-Eyges transport. 

Figures and tables specifically for this paper were computed with various branches of SPWR.FOR. 
ProjScatPower(x) computes the projected scattering power at a; in a stack according to the formula 
selected by software switch spMode. WhatsHere(x...) returns material properties and such quanti- 
ties as Kanematsu's (. and Schneider's t at x. The body of each table set in I^TJ^X was produced as 
a text file to avoid errors of transcription. Data matrices required for graphics were also imported 
from text files. 

We computed isolated examples with our 'proton desk calculator' LOOKUP, a WinXP executable 
distributed with BGware. LOOKUP is a convenient driver for some of the subroutines mentioned, 
offering a choice of 'tasks'. A useful one in the present context is STACK, which computes various 
quantities as a proton beam proceeds down a stack of slabs. 

Results depend somewhat on the choice of range-energy tables. Common tables differ by 1- 
2% for a given material, presumably because of different choices of /, the mean excitation energy. 
We used ICRU Report 49 throughout (our table MIXED. RET in \BGware\data) despite some 
experimental evidence [27| that Janni's 1982 tables [5S] are better for water. Reference [4j used 
Janni's 1966 tables ^5] , 

Finally, some comments on integration. We need to evaluate (Eq.[2]) at many values oi x/Ri 
for the numerous graphs of 9xx compared with ^Hanson- By far the most efficient way is to divide 
the slab by equal ratios so that the contributions of each step to the sum are nearly equal instead of 
very different r^l If, in addition, Simpson's rule is used rather than the midpoint rule, extremely fast 
and accurate integration is achieved. The math (subroutine SimpRat in module SPWRSUBS.FOR) 
is slightly confusing. Suppose we wish to divide a slab of thickness x into n steps (Ax)i . . . (Ax)n 
such that the ratio of successive steps is a constant r. Let us assume some provisional value for n. 
The formula for the sum of a geometric progression 

l + r + r^ + .-. + r"-! = (r"-l)/(r-l) 

leads directly to formulas for the last and first steps 

{Ax)n = ^— 4 a; , (Ax)i = (Aa:)„ 



■^^ In our code T is kinetic energy, E is total energy, R is in g/cm'^ and the longitudinal coordinate is usually z. 
In general, we interpolate published tables for range-energy relations but compute MCS quantities directly from 
one of the applicable theories. We use few analytical approximations, but many lookup tables. 
^® There is no advantage in dividing the higher Fermi-Eyges integrals Ai and A2 by ratios. 
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and we simply divide by r on each iteration. That works for any integer n > 1 and any r > 1. It is 
physicaUy reasonable to let 



Ri 



Ri — X 



l/n 



because the need to subdivide at all (variation of pv) occurs only when the residual range (i?i — x) is 
significantly smaller than the range Ri. By the same argument it is foolish to pick n at the outset, 
since that leads to unnecessary subdivision of thin slabs. Therefore we input some desired maximum 
value of r and let 

INT(l.+ ^^^(^^/(^^-^))~ 



In(rmax) 

For Ri — X = 0.03 i?i and r^ax — 1-6 we find n ~ 8, the most steps needed in practice. 

Unlike the midpoint rule, Simpson's rule uses r(0) where any nonlocal formula for T diverges 
because the single scattering correction diverges. That is easily fixed by substituting a very small 
positive value, a given fraction of the first step size, for 0. Final results are quite insensitive to that 
fraction, which we have adjusted to suppress discontinuities, as n takes on successive values, in the 
graphs of 9xx vs. x/Ri. 



B Preston and Koehler's Formula 

Preston and Koehler [23j show that, to a good approximation, the rms radius a = \pl yrms at any 
normalized depth t in any single slab is related to its maximum value co by 



2(1 -i)2ln 



1 



1-t 



3r 



2t 



1/2 



where t = x/ Ri is the normalized depth, uo in water is 

ao = 0.00627 F^/^ Ri "'^^^ cm 

where 



w^Z,{Z, + l)_^^^ r 106 



A, 



1/2 



(41) 

(42) 
(43) 



and /3 is given by ([5]) with E corresponding to the kinetic energy at a depth Ri/2. 

Their derivation of (j4ip involves the 0veras approximation which they apparently discovered 
independently. They use a quantity similar to Tic based on the multiple scattering formalism of 
Bethe and Ashkin f5D]. Kanematsu [H] gives a more modern derivation of (|^T|) based on Tld- 
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4 Orms according to eight models for various normalized thicknesses and materials. From 
left: normalized slab thickness, actual slab thickness, outgoing energy, ^Hanson in 
milliradians. Remaining seven entries are difference from ^Hanson expressed in %, 
100 X (0xx/feanson — 1)- ^LD is includcd for reference only; it is not supposed to be 
valid for non tissue like matter 

5 (T/p) according to eight models for various materials and normalized thicknesses. 
From left: normalized slab thickness, actual slab thickness, outgoing energy, (r//9)Hansc 
in milliradian^/(g/cm^). Remaining seven entries are difference from (r/p)Hanson ex- 
pressed in %, 100 X {{T/p)xx/ (T/p) Hanson — 1)- {T / p)ld is included for reference 
only; it is not supposed to be valid for non tissue like matter 

6 Pb/Lexan combinations of the simplified range modulator used in computing Figure 
[201 Ei is the proton energy entering the i^^ slab {Ei = 230 MeV), 6'Hanson,3 is the 
design MCS angle entering air and xq is the effective scattering point used in designing 
the modulator for constant yrms = 3.5 cm at 100 cm. The puUback per position is 
2.308 cm water equivalent 
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158.60 MeV p on Be, pRi = 21.290 g/cm^ : 



x/Rx 


px 


E 


^Hanson 


^^Highland 








00S 


^^dH 


^dM 




g/cm^ 


MeV 


mrad 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.021 


158.51 


0.565 


-6.15 


62.86 


36.65 


45.63 


12.10 


-6.51 


-3.03 


0.010 


0.213 


157.69 


2.020 


-1.62 


44.55 


21.29 


29.23 


0.59 


-1.89 


-0.22 


0.100 


2.129 


149.32 


7.207 


3.16 


31.43 


10.28 


17.30 


-1.52 


3.10 


1.63 


0.200 


4.258 


139.62 


10.784 


4.56 


28.12 


7.50 


14.11 


-0.43 


4.72 


2.02 


0.500 


10.645 


107.00 


19.812 


6.31 


23.93 


3.98 


9.45 


-2.17 


7.18 


2.21 


0.900 


19.161 


43.71 


37.994 


7.35 


21.34 


1.82 


3.94 


-8.96 


9.81 


1.14 


0.970 


20.651 


22.49 


48.051 


7.76 


21.33 


1.80 


1.27 


-13.65 


11.00 


0.14 


158.60 MeV p on Al, 


= 22.372 j 


l/cw? : 














x/Ri 


px 


E 


^Hanson 


^^Highland 










^^dH 


^dM 




g/cni^ 


McV 


mracl 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.022 


158.51 


1.004 


-3.33 


54.97 


41.63 


38.57 


14.29 


-3.67 


0.56 


0.010 


0.224 


157.68 


3.658 


-1.82 


34.83 


23.22 


20.54 


1.01 


-2.08 


1.43 


0.100 


2.237 


149.22 


13.214 


0.82 


21.11 


10.69 


8.07 


0.54 


0.77 


2.06 


0.200 


4.474 


139.42 


19.836 


1.68 


17.72 


7.59 


4.80 


2.61 


1.83 


2.14 


0.500 


11.186 


106.52 


36.631 


2.72 


13.45 


3.68 


0.06 


0.93 


3.53 


1.94 


0.900 


20.135 


42.99 


70.843 


3.28 


10.78 


1.25 


-5.39 


-7.92 


5.56 


0.56 


0.970 


21.701 


21.83 


90.068 


3.69 


10.81 


1.27 


-7.50 


-14.43 


6.70 


-0.45 


158.60 MeV p on Cu, pRi 


= 26.258 g/cm^ : 














x/Ri 


px 


E 


^Hanson 


highland 


^FR 


Oic 




00S 


^dH 


^dM 




g/cni^ 


McV 


mrad 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.026 


158.51 


1.545 


-1.71 


49.12 


39.85 


33.34 


11.88 


-1.99 


-0.60 


0.010 


0.263 


157.67 


5.664 


-1.56 


28.90 


20.88 


15.23 


-1.66 


-1.79 


-0.44 


0.100 


2.626 


149.14 


20.535 


0.17 


15.40 


8.22 


2.95 


-1.80 


0.13 


-0.18 


0.200 


5.252 


139.25 


30.855 


0.81 


12.10 


5.13 


-0.24 


0.46 


0.97 


-0.15 


0.500 


13.129 


106.08 


57.059 


1.65 


8.02 


1.31 


-4.83 


-0.95 


2.43 


-0.36 


0.900 


23.632 


42.22 


110.652 


2.44 


5.85 


-0.73 


-9.79 


-9.29 


4.66 


-1.41 


0.970 


25.470 


21.14 


140.987 


3.27 


6.32 


-0.29 


-10.85 


-15.57 


6.20 


-2.07 


158.60 MeV p on Pb, pRi 


= 36.057 g/cm2 : 














x/Ri 


px 


E 


^Hanson 


feighland 


^FR 


^IC 






^dH 


^^dM 




g/cm^ 


MeV 


mrad 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.036 


158.51 


2.638 


2.55 


45.33 


42.50 


29.95 


10.10 


2.24 


1.39 


0.010 


0.361 


157.65 


9.878 


-0.31 


23.07 


20.68 


10.02 


-5.13 


-0.53 


-0.51 


0.100 


3.606 


148.96 


36.222 


-0.35 


8.99 


6.87 


-2.80 


-5.95 


-0.38 


-1.35 


0.200 


7.211 


138.91 


54.568 


-0.10 


5.66 


3.60 


-6.05 


-3.83 


0.06 


-1.53 


0.500 


18.029 


105.27 


101.274 


0.41 


1.71 


-0.27 


-10.62 


-5.13 


1.17 


-1.86 


0.900 


32.452 


40.97 


197.440 


1.56 


0.15 


-1.80 


-15.15 


-12.62 


3.71 


-2.51 


0.970 


34.976 


19.95 


253.086 


2.92 


1.15 


-0.82 


-14.93 


-18.53 


5.78 


-2.74 



Table 4: 9rms according to eight models for various normalized thicknesses and materials. From left: 
normalized slab thickness, actual slab thickness, outgoing energy, ^^Hanson in milliradians. Remaining 
seven entries are difference from ^?Hanson expressed in %, 100 x (f?xx/^^Hanson — 1)- ^ld is included 
for reference only; it is not supposed to be valid for non tissue like matter. 
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158.60 MeV p on Be, pRi = 21.290 g/cm^ : 



x/Ri 


px 


E 


(T'/p)Hans 

mr^cm^/g 


Highland 


FR 


IC 


LD 


0S 


dH 


dM 




g/cw? 


MeV 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.021 


158.51 


16.70 


-9.04 


138.70 


68.05 


90.84 


13.36 


-8.42 


-2.16 


0.010 


0.213 


157.69 


20.99 


0.47 


91.78 


35.02 


53.23 


-5.03 


0.60 


1.94 


0.100 


2.129 


149.32 


27.65 


10.18 


61.23 


13.51 


27.96 


-0.36 


10.81 


4.75 


0.200 


4.258 


139.62 


32.89 


12.84 


53.74 


8.24 


21.00 


-0.38 


13.85 


4.79 


0.500 


10.645 


107.00 


57.19 


15.46 


46.27 


2.98 


11.29 


-8.86 


19.07 


4.08 


0.900 


19.161 


43.71 


324.65 


16.67 


45.37 


2.35 


-2.40 


-36.96 


25.45 


-1.16 


0.970 


20.651 


22.49 


1168.89 


18.80 


49.29 


5.11 


-10.78 


-38.70 


29.78 


-4.89 



158.60 MeV p on Al, pR^ = 22.372 g/cm^ : 



x/Ri 


px 


E 


{T/p)]ians 

mr cm /g 


Highland 


FR 


IC 


LD 


0S 


dH 


dM 




g/cm^ 


McV 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.022 


158.51 


51.08 


-6.29 


111.92 


77.00 


69.43 


15.67 


-5.67 


3.21 


0.010 


0.224 


157.08 


66.29 


-1.87 


64.90 


37.73 


31.75 


-4.44 


-1.73 


4.09 


0.100 


2.237 


149.22 


89.29 


3.83 


35.73 


13.37 


7.63 


6.25 


4.33 


4.69 


0.200 


4.474 


139.42 


106.54 


5.43 


29.21 


7.92 


1.51 


7.56 


6.59 


4.55 


0.500 


11.186 


106.52 


187.74 


6.85 


22.01 


1.91 


-7.64 


-3.46 


10.10 


3.01 


0.900 


20.135 


42.99 


1095.53 


7.68 


20.87 


0.95 


-19.14 


-42.14 


15.29 


-2.65 


0.970 


21.701 


21.83 


4015.91 


10.48 


25.13 


4.51 


-21.49 


-45.01 


20.19 


-5.72 



158.60 MeV p on Cu, pRi = 26.258 g/cm^ : 



x/Ri 


px 


E 


(7'//0)Hans 

mr^cm^/g 


Highland 


FR 


IC 


LD 


0S 


dH 


dM 




g/cw? 


McV 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.026 


158.51 


103.73 


-4.34 


94.86 


71.38 


55.79 


10.10 


-3.73 


0.04 


0.010 


0.263 


157.67 


135.90 


-2.17 


50.20 


32.10 


19.99 


-9.53 


-2.03 


-0.07 


0.100 


2.626 


149.14 


183.84 


1.93 


23.21 


8.37 


-2.36 


2.07 


2.57 


0.14 


0.200 


5.252 


139.25 


219.98 


3.18 


17.11 


3.00 


-8.12 


3.61 


4.34 


-0.16 


0.500 


13.129 


106.08 


389.19 


4.57 


10.76 


-2.58 


-16.49 


-6.53 


7.59 


-1.51 


0.900 


23.632 


42.22 


2292.33 


7.02 


11.71 


-1.75 


-25.08 


-43.41 


14.48 


-5.40 


0.970 


25.470 


21.14 


8491.42 


12.04 


17.78 


3.59 


-19.18 


-45.43 


21.51 


-6.88 



158.60 MeV p on Pb, pRi = 36.057 g/cm^ : 



x/Ri 


px 


E 


{T/ p)}iaas 

mi cm /g 


Highland 


FR 


IC 


LD 


0S 


dH 


dM 




g/cvo? 


McV 


% 


% 


% 


% 


% 


% 


% 


0.001 


0.036 


158.51 


225.07 


1.16 


81.30 


74.31 


44.94 


4.47 


1.78 


1.99 


0.010 


0.361 


157.65 


304.43 


-1.40 


35.39 


30.17 


8.15 


-16.64 


-1.25 


-1.35 


0.100 


3.606 


148.96 


419.99 


-0.19 


9.12 


4.91 


-13.67 


-6.69 


0.30 


-2.93 


0.200 


7.211 


138.91 


505.18 


0.51 


3.43 


-0.56 


-19.12 


-5.32 


1.43 


-3.51 


0.500 


18.029 


105.27 


898.09 


1.81 


-1.67 


-5.46 


-26.56 


-13.88 


4.72 


-4.40 


0.900 


32.452 


40.97 


5392.73 


6.85 


1.68 


-2.24 


-32.37 


-46.40 


13.96 


-6.13 


0.970 


34.976 


19.95 


20673.75 


14.12 


9.48 


5.26 


-14.95 


-48.30 


23.50 


-5.96 



Table 5: (T/p) according to eight models for various materials and normalized thicknesses. 
From left: normalized slab thickness, actual slab thickness, outgoing energy, (T/p)Hanson in 
milliradian^/(g/cm^). Remaining seven entries are difference from (T/(o)Hanson expressed in %, 
100 X ((T/p) XX /(T'/p) Hanson — !)• {T / p)hD is included for reference only; it is not supposed to be 
valid for non tissue like matter. 



21 



List of Figures 

1 Mixed slab (stack) geometry. 



2 Weak 0veras approximation: the function 1 — ^ as a function of normalized 
slab thickness for protons of energies 20, 50, 100 and 200 MeV exiting single slabs of 
Be, Cu and Pb (all superimposed). 20 MeV protons leaving a very thin Pb slab 
{x/Ri ~ 0.001) have the largest deviation from the ideal, 20% [2^ 

3 Projected mass scattering powers T/p at 20Mev in Be, Cu and Pb vs. normalized 
overlying slab thickness x/Ri. Formulas for Tfr, Tic and Tnanson are given in the text. [24 

4 Deviation from ^Hanson of ^Highland Computed from the generalized Highland formula 
(|13p . for four scattering materials. The incident proton energy is 158.6 MeV [1^ 

5 Deviation from ^Hanson of 0fr computed from Tfr, the Fermi-Rossi scattering power, 
for four scattering materials. The incident proton energy is 158.6 MeV 

6 Deviation from ^Hanson of 9ic computed from Tic, the ICRU35 Report 35- scattering 
power adapted to protons, for four scattering materials. The incident proton energy 

is 158.6 MeV and the line width code the same as Figure[5l [l^ 

7 Deviation from ^Hanson of Old computed from Tld, Kanematsu's 'linear displacement' 
scattering power [21. , for four scattering materials. The incident proton energy is 
158.6 MeV M 

8 Deviation from ^Hanson of 00g computed from T0s, the scattering power of Schneider 

et al. [TFi, for four scattering materials. The incident proton energy is 158.6 MeV. . |27 

9 Deviation from ^Hanson of 0dH computed from Tdfj, Kanematsu's 'differential High- 
land' scattering power [7], for four scattering materials. The incident proton energy 

is 158.6 MeV [27 

10 Ratio THanson/21c as a function of 1 — (pv/pivi)'^ at four proton energies for three 
materials: Be (small circles), Cu (medium circles) and Pb (large circles). The heavy 
lines are a bilinear fit to the entire data set [l^ 

11 Deviation from ^^Hanson of O^m computed from TdM, the scattering power proposed 
in the present work, for four scattering materials. The incident proton energy is 
158.6 MeV M 

12 For polystyrene, deviation of 0xx from ^Hanson at 158.6 Mev incident energy. Data 
(open circles) taken from 4 , . Each 0xx is the integral of the corresponding Txx except 
^'Highland which is from the generalized Highland formula p^ . 6'fr is off scale. ... 28 

13 The same as Figure [T2| for Pb [26 

14 Spreading of a 127 MeV proton pencil beam in water with experimental data from 
Preston and Koehler [23]. Tho is taken from Table 2 of HoUmark et al. [6]. Apart 
from Tfr and ThOi calculations based on the other scattering powers are barely 
distinguishable [s^ 

15 Expanded view of Figure [14] using Preston and Koehler's analytic approximation to 
beam spreading (see text) as a an arbitrary reference. Line styles are Tfr light solid. 
Tic dotted, Tld dot dash, Tqs short dash, TdH long dash, TdM bold solid for this and 

all following figures 3C 

16 Integrand of A2 a,t x = 0.97 i?i for 127 Mev protons incident on water 31 

17 Same as Figure [THI except x = 0.1 -Ri [31 

18 Beam phase space ellipses for 127 Mev protons incident on water using Fermi-Eyges 
theory with six scattering powers, at the exit faces of five equal slabs extending to 
0.97 i?i [32 

19 Same as Figure [11] for Pb E 

20 Projected rms displacement yxx, according to Fermi-Eyges theory with six formulas 
for T, at the end of a sequence of Pb/Lexan/air stacks corresponding to a simplified 
range modulator (Table[6]). The horizontal line at yxx — 3.5 cm represents the design 

goal [13 



22 



;;R;;; 



Figure 1: Mixed slab (stack) geometry. 
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Figure 2: Weak 0vcras approximation: the function 1 — (pv)'^ /{pv)i ^ as a function of normalized 
slab thickness for protons of energies 20, 50, 100 and 200 MeV exiting single slabs of Be, Cu and Pb 
(all superimposed). 20 MeV protons leaving a very thin Pb slab {x/R\ = 0.001) have the largest 
deviation from the ideal, 20%. 
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Figure 3: Projected mass scattering powers T/ p at 20 Mev in Be, Cu and Pb vs. normalized overlying 
slab thickness x/Ri. Formulas for Tfr, Tic and Tnanson are given in the text. 
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Figure 5: Deviation from ^Hanson of flpR computed from Tfr , the Fermi-Rossi scattering power, for 
four scattering materials. The incident proton energy is 158.6 MeV. 
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Figure 6: Deviation from 6'Hanson of 9ic computed from Tic, the ICRU35 Report 35- scattering 
power adapted to protons, for four scattering materials. The incident proton energy is 158.6 MeV 
and the hne width code the same as Figure[Sl 
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Figure 7: Deviation from ^Hanson of 9ld computed from Tld, Kanematsu's 'hnear displacement' 
scattering power [3T], for four scattering materials. The incident proton energy is 158.6 MeV. 
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Figure 8: Deviation from 6'Hanson of Oqs computed from T0s, the scattering power of Schneider et 
al. [15], for four scattering materials. The incident proton energy is 158.6 MeV. 
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Figure 9: Deviation from ^Hanson of Odn computed from TdH, Kanematsu's 'differential Highland' 
scattering power 7J, for four scattering materials. The incident proton energy is 158.6 MeV. 
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Figure 11: Deviation from ^Hanson of ^dM computed from T^m, the scattering power proposed in the 
present work, for four scattering materials. The incident proton energy is 158.6 MeV. 



28 




0.001 0.01 0.1 1 

X/ 

Figure 12: For polystyrene, deviation of 0xx from ^Hanson at 158.6 Mev incident energy. Data (open 
circles) taken from [3]. Each 0xx is the integral of the corresponding Txx except 6'Highiand which is 
from the generalized Highland formula (fT^l) . 6'fr is off scale. 




Figure 13: The same as Figure [T^ for Pb. 
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Figure 14: Spreading of a 127 MeV proton pencil beam in water with experimental data from 
Preston and Koehler [23| . Tho is taken from Table 2 of HoUmark et al. ,6] . Apart from TpR and 
Tho, calculations based on the other scattering powers are barely distinguishable. 
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Figure 15: Expanded view of Figure [141 using Preston and Koehler's analytic approximation to beam 
spreading (see text) as a an arbitrary reference. Line styles are TpR light solid, Tic dotted, Tld dot 
dash, Tqq, short dash, TdH long dash, T^m bold solid for this and all following figures. 
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Figure 16: Integrand of A2 at x = 0.97 i?i for 127 Mev protons incident on water. 
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Figure 17: Same as Figure [TBI except x — 0.1 Ri. 



31 



80 




-3-2-10 1 2 3 
y (mm) 



Figure 18: Beam phase space ellipses for 127 Mev protons incident on water using Fermi-Eyges 
theory with six scattering powers, at the exit faces of five equal slabs extending to 0.97 Ri. 




Figure 19: Same as Figure [TSl for Pb. 
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Pb / Lexan combination 



Figure 20: Projected rms displacement yxx, according to Fermi-Eyges theory with six formulas for 
T, at the end of a sequence of Pb/Lexan/air stacks corresponding to a simplified range modulator 
(TableO. The horizontal line at yxx = 3.5cm represents the design goal. 



comb 


Pb 


Lexan 


E2 


# 


g/cm^ 


g/cm^ 


MeV 


1 


6.429 


0.000 


216.4 


2 


6.173 


2.560 


216.9 


3 


5.872 


5.144 


217.6 


4 


5.543 


7.743 


218.3 


5 


5.179 


10.360 


219.1 


6 


4.781 


12.995 


219.9 


7 


4.335 


15.656 


220.9 


8 


3.834 


18.346 


221.9 


9 


3.240 


21.085 


223.2 


10 


2.509 


23.898 


224.7 


11 


1.537 


26.840 


226.8 


12 


0.000 


30.082 


72.2 



E3 


E4 


^Hanson, 3 


Xq 


MeV 


MeV 


mrad 


cm 


216.4 


215.9 


35.10 


0.28 


206.4 


205.9 


35.12 


0.34 


196.1 


195.6 


35.19 


0.53 


185.4 


184.9 


35.32 


0.88 


174.3 


173.8 


35.51 


1.44 


162.7 


162.2 


35.79 


2.23 


150.5 


150.0 


36.19 


3.30 


137.5 


137.0 


36.73 


4.71 


123.7 


123.1 


37.45 


6.55 


108.6 


108.0 


38.42 


8.94 


91.8 


91.1 


39.78 


12.11 


72.2 


71.5 


41.93 


16.54 



Table 6: Pb/Lexan combinations of the simplified range modulator used in computing Figure [20l 
Ei is the proton energy entering the j"^ slab {Ei = 230 MeV), 6'Hanson, 3 is the design MCS angle 
entering air and xq is the effective scattering point used in designing the modulator for constant 
Vrms — 3.5 cm at 100 cm. The pullback per position is 2.308 cm water equivalent. 
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